{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "inputHidden": false,
    "outputHidden": false
   },
   "outputs": [],
   "source": [
    "from IPython.core.interactiveshell import InteractiveShell\n",
    "InteractiveShell.ast_node_interactivity = \"all\"\n",
    "from functools import *\n",
    "\n",
    "import sys\n",
    "sys.path.append('../GTF')\n",
    "\n",
    "import networkx as nx\n",
    "import pygsp as pg\n",
    "import numpy as np\n",
    "import random\n",
    "import scipy as sp\n",
    "\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "from IPython.display import display\n",
    "from math import log\n",
    "#from tqdm import tqdm, trange\n",
    "from tqdm import tqdm_notebook\n",
    "from tqdm import tnrange\n",
    "import math"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "inputHidden": false,
    "outputHidden": false
   },
   "outputs": [],
   "source": [
    "#create container of num_iter piecewise graph signals:\n",
    "def noisy_pwc(num_iter, Gnx, num_seeds, noise_level):\n",
    "    beta_0_list = []\n",
    "    y_list = []\n",
    "    path_lens = dict(nx.shortest_path_length(Gnx))\n",
    "    for it in range(num_iter):\n",
    "        beta_0 = gen_pwc_sig(Gnx,num_seeds,path_lens)\n",
    "        y = beta_0.copy() + noise_level * np.random.rand(Gnx.number_of_nodes())\n",
    "        beta_0_list.append(beta_0)\n",
    "        y_list.append(y)\n",
    "    \n",
    "    beta_0_list = np.array(beta_0_list)\n",
    "    y_list = np.array(y_list)\n",
    "    return {'beta_0':beta_0_list,'y':y_list}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "inputHidden": false,
    "outputHidden": false
   },
   "outputs": [],
   "source": [
    "def controller(gamma, Gnx, k, rho, penalty_param,penalty_f,signals, num_iter, beta_init):\n",
    "    mse_rs = 0.\n",
    "#     l0diff_rs = 0.\n",
    "    l0beta_rs = 0.\n",
    "    suppTPR_rs = 0.\n",
    "    suppFPR_rs = 0.\n",
    "    betas = []\n",
    "    for it in range(num_iter):\n",
    "        beta_0 = signals['beta_0'][it]\n",
    "        y = signals['y'][it]      \n",
    "        if beta_init is not None:\n",
    "            out = admm_mse2(gamma,rho,penalty_param,Gnx,beta_0,y,k,penalty_f,\n",
    "                           beta_init=beta_init[it])\n",
    "        else:\n",
    "            out = admm_mse2(gamma,rho,penalty_param,Gnx,beta_0,y,k,penalty_f,\n",
    "                           beta_init=None)          \n",
    "        mse_rs += out['mse']\n",
    "#         l0diff_rs += out['l0diff']\n",
    "        l0beta_rs += out['l0beta']\n",
    "        suppTPR_rs += out['suppTPR']\n",
    "        suppFPR_rs += out['suppFPR']\n",
    "        betas.append(out['beta'])\n",
    "    \n",
    "    output = {'avg_mse':mse_rs/num_iter, #'avg_l0diff':np.round(l0diff_rs/num_iter), \n",
    "            'avg_l0beta':np.around(l0beta_rs/num_iter, decimals=3), 'avg_suppTPR':np.around(suppTPR_rs/num_iter, decimals=3),\n",
    "           'avg_suppFPR':np.around(suppFPR_rs/num_iter, decimals=3)}  \n",
    "    if penalty_f == 'L1':\n",
    "        output['beta'] = betas\n",
    "\n",
    "    return output"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "inputHidden": false,
    "outputHidden": false
   },
   "outputs": [],
   "source": [
    "G = pg.graphs.Minnesota(connect = True)\n",
    "#G = pg.graphs.Grid2d(20,20)\n",
    "#G = pg.graphs.ErdosRenyi(100,0.05)\n",
    "#G = pg.graphs.Sensor(20,1)\n",
    "#G.set_coordinates()\n",
    "Gnx = nx.from_scipy_sparse_matrix(G.A.astype(float),edge_attribute = 'weight')\n",
    "num_seeds = 4\n",
    "path_lens = dict(nx.shortest_path_length(Gnx))\n",
    "beta_0 = gen_pwc_sig(Gnx,num_seeds, path_lens) \n",
    "#beta_0 = beta_0 / np.linalg.norm(beta_0,2)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "LOW NOISE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "inputHidden": false,
    "outputHidden": false
   },
   "outputs": [],
   "source": [
    "noise_level = 0.05\n",
    "\n",
    "y = beta_0 + noise_level * np.random.rand(G.N)\n",
    "limits = [ min(y),  max(y)]\n",
    "#example beta_0,y for current parameter setting\n",
    "\n",
    "fig, axes = plt.subplots(nrows=1, ncols=2,figsize=(15, 6))\n",
    "pg.plotting.plot_signal(G,beta_0,limits = limits, plot_name = \"Beta_0\", ax = axes[0])\n",
    "pg.plotting.plot_signal(G,y,limits = limits,plot_name = \"y : noisy\", ax = axes[1])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "inputHidden": false,
    "outputHidden": false
   },
   "outputs": [],
   "source": [
    "num_iter = 1\n",
    "signals = noisy_pwc(num_iter, Gnx, num_seeds, noise_level )\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "inputHidden": false,
    "outputHidden": false
   },
   "outputs": [],
   "source": [
    "\n",
    "num_pts = 100\n",
    "log_gamma_sweep = np.linspace(-3,3,num_pts)\n",
    "gamma_sweep = [10**(y) for y in log_gamma_sweep]\n",
    "\n",
    "num_pts2 = 20\n",
    "log_rhomult_sweep = np.linspace(-1,1,num_pts2)\n",
    "rhomult_sweep = [10**(y) for y in log_rhomult_sweep]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "inputHidden": false,
    "outputHidden": false
   },
   "outputs": [],
   "source": [
    "#find best for given gamma by sweeping over rho\n",
    "import copy \n",
    "def sweep_rho(rm_sweep, gamma,Gnx,k,penalty_param,num_iter,penalty_f, signals , beta_init):\n",
    "    min_mse = float('inf')\n",
    "    \n",
    "    for jj in tnrange(len(rm_sweep)):   \n",
    "        #print jj\n",
    "        wrapper_gtf = controller(gamma, Gnx=Gnx,k=k,rho=max(rm_sweep[jj]*gamma,1.0/penalty_param),penalty_param=penalty_param,\n",
    "                        num_iter=num_iter,penalty_f = penalty_f,signals = signals, beta_init=beta_init)\n",
    "        mse = wrapper_gtf['avg_mse']\n",
    "        if(mse < min_mse):\n",
    "            min_mse = mse \n",
    "            wrapper_min = copy.deepcopy(wrapper_gtf)    \n",
    "    return wrapper_min\n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "inputHidden": false,
    "outputHidden": false
   },
   "outputs": [],
   "source": [
    "k = 0\n",
    "out_L1 = []\n",
    "out_SCAD = []\n",
    "out_MCP = []\n",
    "for ii in  tnrange(num_pts):\n",
    "    print ii\n",
    "    wrapper_L1 = sweep_rho(rm_sweep = rhomult_sweep, gamma = gamma_sweep[ii],Gnx = Gnx,k = k\n",
    "                           ,penalty_param = float('inf') ,num_iter = num_iter,penalty_f = \"L1\", signals = signals , beta_init = None)\n",
    "    wrapper_SCAD = sweep_rho(rm_sweep = rhomult_sweep, gamma = gamma_sweep[ii],Gnx = Gnx,k = k\n",
    "                           ,penalty_param = 3.7 ,num_iter = num_iter,penalty_f = \"SCAD\", signals = signals , beta_init = wrapper_L1['beta'])\n",
    "    wrapper_MCP = sweep_rho(rm_sweep = rhomult_sweep, gamma = gamma_sweep[ii],Gnx = Gnx,k = k\n",
    "                           ,penalty_param = 1.4 ,num_iter = num_iter,penalty_f = \"MCP\", signals = signals , beta_init = wrapper_L1['beta'])\n",
    "\n",
    "    del wrapper_L1['beta']\n",
    "    out_L1.append(wrapper_L1)\n",
    "    out_SCAD.append(wrapper_SCAD)\n",
    "    out_MCP.append(wrapper_MCP)\n",
    "        \n",
    "        \n",
    "        \n",
    "        \n",
    "out_L1 = {k: [dic[k] for dic in out_L1] for k in out_L1[0] if k is not 'beta'}\n",
    "out_SCAD = {k: [dic[k] for dic in out_SCAD] for k in out_SCAD[0]}\n",
    "out_MCP = {k: [dic[k] for dic in out_MCP] for k in out_MCP[0]}\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "inputHidden": false,
    "outputHidden": false
   },
   "outputs": [],
   "source": [
    "ind = np.argsort(out_L1['avg_l0beta'])\n",
    "x_L1 = np.array([out_L1['avg_l0beta'][i] for i in ind]).astype(int)\n",
    "mse_L1 = np.array([out_L1['avg_mse'][i] for i in ind])\n",
    "# l0diff_L1 = np.array([out_L1['avg_l0diff'][i] for i in ind])\n",
    "ind2 = np.argsort(out_L1['avg_suppFPR'])\n",
    "suppTPR_L1 = np.array([out_L1['avg_suppTPR'][i] for i in ind2])\n",
    "suppFPR_L1 = np.array([out_L1['avg_suppFPR'][i] for i in ind2])\n",
    "\n",
    "ind = np.argsort(out_SCAD['avg_l0beta'])\n",
    "x_SCAD = np.array([out_SCAD['avg_l0beta'][i] for i in ind]).astype(int)\n",
    "mse_SCAD = np.array([out_SCAD['avg_mse'][i] for i in ind])\n",
    "# l0diff_SCAD = np.array([out_SCAD['avg_l0diff'][i] for i in ind])\n",
    "ind2 = np.argsort(out_SCAD['avg_suppFPR'])\n",
    "suppTPR_SCAD = np.array([out_SCAD['avg_suppTPR'][i] for i in ind2])\n",
    "suppFPR_SCAD = np.array([out_SCAD['avg_suppFPR'][i] for i in ind2])\n",
    "\n",
    "ind = np.argsort(out_MCP['avg_l0beta'])\n",
    "x_MCP = np.array([out_MCP['avg_l0beta'][i] for i in ind]).astype(int)\n",
    "mse_MCP = np.array([out_MCP['avg_mse'][i] for i in ind])\n",
    "# l0diff_MCP = np.array([out_MCP['avg_l0diff'][i] for i in ind])\n",
    "ind2 = np.argsort(out_MCP['avg_suppFPR'])\n",
    "suppTPR_MCP = np.array([out_MCP['avg_suppTPR'][i] for i in ind2])\n",
    "suppFPR_MCP = np.array([out_MCP['avg_suppFPR'][i] for i in ind2])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "inputHidden": false,
    "outputHidden": false
   },
   "outputs": [],
   "source": [
    "#mse plot:\n",
    "plt.figure()\n",
    "plt.plot(x_L1, mse_L1, label='L1')\n",
    "plt.plot(x_SCAD, mse_SCAD, label='SCAD')\n",
    "plt.plot(x_MCP, mse_MCP, label='MCP')\n",
    "plt.legend()\n",
    "\n",
    "plt.figure()\n",
    "plt.plot(suppFPR_L1, suppTPR_L1, 'o--', label='L1')\n",
    "plt.plot(suppFPR_SCAD, suppTPR_SCAD, 'v--', label='SCAD')\n",
    "plt.plot(suppFPR_MCP, suppTPR_MCP, '.--', label='MCP')\n",
    "plt.xlabel('False Positive Rate')\n",
    "plt.ylabel('True Positive Rate')\n",
    "plt.legend()\n",
    "plt.title('Support Recovery ROC Curve')\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "inputHidden": false,
    "outputHidden": false
   },
   "outputs": [],
   "source": [
    "def smooth_curve(bins, x, y):\n",
    "    mask = np.digitize(x, bins) # bins[i-1]< x <= bins[i]\n",
    "    return np.array([y[mask==idx].mean() for idx in range(1,len(bins))])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "inputHidden": false,
    "outputHidden": false
   },
   "outputs": [],
   "source": [
    "#binning (plot average value in each bin) to smooth out the curve\n",
    "#print signals['y'].size\n",
    "#print len(Gnx.edges)\n",
    "#print x_L1\n",
    "#print x_SCAD\n",
    "#print x_MCP\n",
    "num_bins = 20\n",
    "bins_mse = np.linspace(min(x_L1.min(), x_SCAD.min(), x_MCP.min()), max(x_L1.max(), x_SCAD.max(), x_MCP.max()), num_bins)\n",
    "#print bins\n",
    "mse_L1_avg = smooth_curve(bins_mse, x_L1, mse_L1)\n",
    "mse_SCAD_avg = smooth_curve(bins_mse, x_SCAD, mse_SCAD)\n",
    "mse_MCP_avg = smooth_curve(bins_mse, x_MCP, mse_MCP)\n",
    "\n",
    "plt.figure()\n",
    "plt.plot(bins_mse[:-1], mse_L1_avg, 'o--', label='L1')\n",
    "plt.plot(bins_mse[:-1], mse_SCAD_avg, 'v--', label='SCAD')\n",
    "plt.plot(bins_mse[:-1], mse_MCP_avg, '.--',label='MCP')\n",
    "plt.legend()\n",
    "plt.title('mse')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "inputHidden": false,
    "outputHidden": false
   },
   "outputs": [],
   "source": [
    "num_bins = 20\n",
    "bins_fpr = np.linspace(min(suppFPR_L1.min(), suppFPR_SCAD.min(), suppFPR_MCP.min()), max(suppFPR_L1.max(), suppFPR_SCAD.max(), suppFPR_MCP.max()), num_bins)\n",
    "#print bins\n",
    "suppTPR_L1_s = smooth_curve(bins_fpr, suppFPR_L1, suppTPR_L1)\n",
    "suppTPR_SCAD_s = smooth_curve(bins_fpr, suppFPR_SCAD, suppTPR_SCAD)\n",
    "suppTPR_MCP_s = smooth_curve(bins_fpr, suppFPR_MCP, suppTPR_MCP)\n",
    "\n",
    "\n",
    "\n",
    "plt.figure()\n",
    "plt.plot(bins_fpr[:-1], suppTPR_L1_s, 'o--', label='L1')\n",
    "plt.plot(bins_fpr[:-1], suppTPR_SCAD_s, 'v--', label='SCAD')\n",
    "plt.plot(bins_fpr[:-1], suppTPR_MCP_s, '.--', label='MCP')\n",
    "plt.xlabel('False Positive Rate')\n",
    "plt.ylabel('True Positive Rate')\n",
    "plt.legend()\n",
    "plt.title('Support Recovery ROC Curve')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "inputHidden": false,
    "outputHidden": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "inputHidden": false,
    "outputHidden": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "inputHidden": false,
    "outputHidden": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "inputHidden": false,
    "outputHidden": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "inputHidden": false,
    "outputHidden": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "inputHidden": false,
    "outputHidden": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "inputHidden": false,
    "outputHidden": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "inputHidden": false,
    "outputHidden": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "inputHidden": false,
    "outputHidden": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "inputHidden": false,
    "outputHidden": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "inputHidden": false,
    "outputHidden": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "inputHidden": false,
    "outputHidden": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "inputHidden": false,
    "outputHidden": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "inputHidden": false,
    "outputHidden": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "HIGH NOISE\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "inputHidden": false,
    "outputHidden": false
   },
   "outputs": [],
   "source": [
    "noise_level= 1\n",
    "y = beta_0 + noise_level * np.random.rand(G.N)\n",
    "limits = [ min(y),  max(y)]\n",
    "\n",
    "\n",
    "#example beta_0,y for current parameter setting\n",
    "fig, axes = plt.subplots(nrows=1, ncols=2,figsize=(15, 6))\n",
    "pg.plotting.plot_signal(G,beta_0,limits = limits, plot_name = \"Beta_0\", ax = axes[0])\n",
    "pg.plotting.plot_signal(G,y,limits = limits,plot_name = \"y : noisy\", ax = axes[1])\n",
    "\n",
    "num_iter = 1\n",
    "num_pts = 1000\n",
    "log_gamma_sweep = np.linspace(-4,8,num_pts)\n",
    "log_gamma_sweep = [10**(y) for y in log_gamma_sweep]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "inputHidden": false,
    "outputHidden": false
   },
   "outputs": [],
   "source": [
    "k = 0\n",
    "out_L1 = []\n",
    "out_SCAD = []\n",
    "out_MCP = []\n",
    "for ii in  tnrange(num_pts):\n",
    "    print ii\n",
    "    wrapper_L1 = controller(log_gamma_sweep[ii], Gnx=Gnx,k=k,rho=log_gamma_sweep[ii],penalty_param=0,#penalty_param,\n",
    "                    num_iter=num_iter,penalty_f = \"L1\",signals = signals, beta_init=None)\n",
    "    wrapper_SCAD = controller(log_gamma_sweep[ii], Gnx=Gnx,k=k,rho=max(log_gamma_sweep[ii], 1.0/3.7),penalty_param=3.7,#penalty_param,\n",
    "                    num_iter=num_iter,penalty_f = \"SCAD\",signals = signals, beta_init=wrapper_L1['beta'])\n",
    "    wrapper_MCP = controller(log_gamma_sweep[ii], Gnx=Gnx,k=k,rho=max(log_gamma_sweep[ii],1.0/1.4),penalty_param=1.4,#penalty_param,\n",
    "                    num_iter=num_iter,penalty_f = \"MCP\",signals = signals, beta_init=wrapper_L1['beta'])\n",
    "    del wrapper_L1['beta']\n",
    "    out_L1.append(wrapper_L1)\n",
    "    out_SCAD.append(wrapper_SCAD)\n",
    "    out_MCP.append(wrapper_MCP)\n",
    "out_L1 = {k: [dic[k] for dic in out_L1] for k in out_L1[0] if k is not 'beta'}\n",
    "out_SCAD = {k: [dic[k] for dic in out_SCAD] for k in out_SCAD[0]}\n",
    "out_MCP = {k: [dic[k] for dic in out_MCP] for k in out_MCP[0]}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "inputHidden": false,
    "outputHidden": false
   },
   "outputs": [],
   "source": [
    "ind = np.argsort(out_L1['avg_l0beta'])\n",
    "x_L1 = np.array([out_L1['avg_l0beta'][i] for i in ind]).astype(int)\n",
    "mse_L1 = np.array([out_L1['avg_mse'][i] for i in ind])\n",
    "# l0diff_L1 = np.array([out_L1['avg_l0diff'][i] for i in ind])\n",
    "ind2 = np.argsort(out_L1['avg_suppFPR'])\n",
    "suppTPR_L1 = np.array([out_L1['avg_suppTPR'][i] for i in ind2])\n",
    "suppFPR_L1 = np.array([out_L1['avg_suppFPR'][i] for i in ind2])\n",
    "\n",
    "ind = np.argsort(out_SCAD['avg_l0beta'])\n",
    "x_SCAD = np.array([out_SCAD['avg_l0beta'][i] for i in ind]).astype(int)\n",
    "mse_SCAD = np.array([out_SCAD['avg_mse'][i] for i in ind])\n",
    "# l0diff_SCAD = np.array([out_SCAD['avg_l0diff'][i] for i in ind])\n",
    "ind2 = np.argsort(out_SCAD['avg_suppFPR'])\n",
    "suppTPR_SCAD = np.array([out_SCAD['avg_suppTPR'][i] for i in ind2])\n",
    "suppFPR_SCAD = np.array([out_SCAD['avg_suppFPR'][i] for i in ind2])\n",
    "\n",
    "ind = np.argsort(out_MCP['avg_l0beta'])\n",
    "x_MCP = np.array([out_MCP['avg_l0beta'][i] for i in ind]).astype(int)\n",
    "mse_MCP = np.array([out_MCP['avg_mse'][i] for i in ind])\n",
    "# l0diff_MCP = np.array([out_MCP['avg_l0diff'][i] for i in ind])\n",
    "ind2 = np.argsort(out_MCP['avg_suppFPR'])\n",
    "suppTPR_MCP = np.array([out_MCP['avg_suppTPR'][i] for i in ind2])\n",
    "suppFPR_MCP = np.array([out_MCP['avg_suppFPR'][i] for i in ind2])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "inputHidden": false,
    "outputHidden": false
   },
   "outputs": [],
   "source": [
    "#mse plot:\n",
    "plt.figure()\n",
    "plt.plot(x_L1, mse_L1, label='L1')\n",
    "plt.plot(x_SCAD, mse_SCAD, label='SCAD')\n",
    "plt.plot(x_MCP, mse_MCP, label='MCP')\n",
    "plt.legend()\n",
    "\n",
    "plt.figure()\n",
    "plt.plot(suppFPR_L1, suppTPR_L1, 'o--', label='L1')\n",
    "plt.plot(suppFPR_SCAD, suppTPR_SCAD, 'v--', label='SCAD')\n",
    "plt.plot(suppFPR_MCP, suppTPR_MCP, '.--', label='MCP')\n",
    "plt.xlabel('False Positive Rate')\n",
    "plt.ylabel('True Positive Rate')\n",
    "plt.legend()\n",
    "plt.title('Support Recovery ROC Curve')\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "inputHidden": false,
    "outputHidden": false
   },
   "outputs": [],
   "source": [
    "#binning (plot average value in each bin) to smooth out the curve\n",
    "\n",
    "bins_mse = np.linspace(min(x_L1.min(), x_SCAD.min(), x_MCP.min()), max(x_L1.max(), x_SCAD.max(), x_MCP.max()), 20)\n",
    "#print bins\n",
    "mse_L1_avg = smooth_curve(bins_mse, x_L1, mse_L1)\n",
    "mse_SCAD_avg = smooth_curve(bins_mse, x_SCAD, mse_SCAD)\n",
    "mse_MCP_avg = smooth_curve(bins_mse, x_MCP, mse_MCP)\n",
    "\n",
    "plt.figure()\n",
    "plt.plot(bins_mse[:-1], mse_L1_avg, 'o--', label='L1')\n",
    "plt.plot(bins_mse[:-1], mse_SCAD_avg, 'v--', label='SCAD')\n",
    "plt.plot(bins_mse[:-1], mse_MCP_avg, '.--',label='MCP')\n",
    "plt.legend()\n",
    "plt.title('mse')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "inputHidden": false,
    "outputHidden": false
   },
   "outputs": [],
   "source": [
    "num_bins = 20\n",
    "bins_fpr = np.linspace(min(suppFPR_L1.min(), suppFPR_SCAD.min(), suppFPR_MCP.min()), max(suppFPR_L1.max(), suppFPR_SCAD.max(), suppFPR_MCP.max()), num_bins)\n",
    "#print bins\n",
    "suppTPR_L1_s = smooth_curve(bins_fpr, suppFPR_L1, suppTPR_L1)\n",
    "suppTPR_SCAD_s = smooth_curve(bins_fpr, suppFPR_SCAD, suppTPR_SCAD)\n",
    "suppTPR_MCP_s = smooth_curve(bins_fpr, suppFPR_MCP, suppTPR_MCP)\n",
    "\n",
    "\n",
    "\n",
    "plt.figure()\n",
    "plt.plot(bins_fpr[:-1], suppTPR_L1_s, 'o--', label='L1')\n",
    "plt.plot(bins_fpr[:-1], suppTPR_SCAD_s, 'v--', label='SCAD')\n",
    "plt.plot(bins_fpr[:-1], suppTPR_MCP_s, '.--', label='MCP')\n",
    "plt.xlabel('False Positive Rate')\n",
    "plt.ylabel('True Positive Rate')\n",
    "plt.legend()\n",
    "plt.title('Support Recovery ROC Curve')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "inputHidden": false,
    "outputHidden": false
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernel_info": {
   "name": "python2"
  },
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.15"
  },
  "nteract": {
   "version": "0.11.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
